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Abstract 



In this paper we study a Hamiltonian system with a spatially asym- 
metric potential. We are interested in the effects on the dynamics when 
the potential becomes symmetric slowly in time. We focus on a highly 
simplified non-trivial model problem (a metaphor) to be able to pur- 
sue explicit calculations as far as possible. Using the techniques of 
averaging and adiabatic invariants, we are able to study all bounded 
solutions, which reveals significant asymmetric dynamics even when 
the asymmetric contributions to the potential have become negligibly 
small. 

1 Introduction 

Many physical objects exhibit some form of symmetry. Most galaxies for 
instance, have axes or planes of symmetry. The motivation for this study is 
that a symmetric equilibrium configuration generally is the outcome of the 
evolution from an asymmetric state. We would like to trace the effect of the 
asymmetries. 

A problem is that studies of the evolution of actual physical systems 
are difficult and so relatively rare. We propose therefore to ignore, at least 
for the time being, the actual physical mechanisms and to consider systems 
described by a Hamiltonian of the form 



where Tis is the part of the Hamiltonian which is symmetric in some 
sense; Tia is the asymmetric part which is slowly vanishing as we put 
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Figure 1: The dynamics of the unperturbed system (^) 



a(0) = 1, lim a{et) = 0, < e < 1 (2) 



To study the dynamics induced by the Hamiltonian Ti.{p, q, et) is still a 
formidable problem. So we simplify as much as possible to obtain 



Xi = X2 

r.2 



I ±2 = —xi + a{et)xi 
which is derived from the one degree of freedom Hamiltonian 



Xi = X2 

±2 = -Xi + xj 



(3) 



nip,q,et) = ^ip^ + q^) + ^ai€t)q^ (4) 

identifying p = X2, q = xi. We shall associate with system the 
"unperturbed" system which arises for e = 



(5) 



We note that in the autonomous system (|5|) there are basically two re- 
gions (figure Q): within the homoclinic solution the orbits are bounded, 
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outside the homochnic solution the orbits diverge to infinity (with the ex- 
ception of the stable manifold and the saddle point itself). In system (^) 
we have no fixed saddle point, still it turns out that we have two separate 
regions of initial values in which the orbits are bounded or diverge to infinity. 

Since the dynamics of systems (^) and (^) are the same on an 0(1) 
timescale, it is instructive (though slightly wrong) to view system (|3|) as 
having a saddle point moving slowly towards infinity and having a slowly ex- 
panding homoclinic orbit. Within this picture, an orbit can remain bounded 
in two ways, either by starting inside the homoclinic orbit, or by getting 
"captured" by the slowly expanding homoclinic orbit, which can only hap- 
pen if the orbit starts sufficiently close to the stable manifold of the saddle 
point. 

Using a special transformation we shall discuss the boundaries of these 
regions in section ^. A special case, a{et) = exp(— et) can be studied easily 
and help us to understand the general case. 

In section |3| we perform averaging in the so-called stable region where 
bounded solutions are found. This involves the use of elliptic and hyperge- 
ometric functions, rather hard analysis, where we are supported by Mathe- 
matica 2.2 running under SunOS 4.1.3. 

After determining the validity of the averaged equation we establish the ex- 
istence of an adiabatic invariant in the stable region, valid for all time. Even 
more remarkable is that explicit calculations of this invariant show that the 
evolution of phase points will show significant traces of its asymmetrical 
past for all time. 

In section ^ we need subtle reasoning to discuss what is going on in the 
boundary layer near the boundary of the stable domain. 

The analysis in this paper is based on averaging methods but, because of 
its direct relation to dissipative mechanics (section ^) , it clearly profits from 
the results by Haberman and Ho |^ Q and Bourland and Haberman |^] . At 
the same time our analysis should be placed within the theory of adiabatic 
invariants, which has been summarized recently in an admirable survey by 
Henrard 

We finally note that in the context of galactic dynamics, some rather 
different examples based on classical results of the theory of adiabatic in- 
variants were given by Binney and Tremaine |l|. 

2 The boundary of the stable part of phase space 

As we explained in the introduction, the phase space of system (P) can 
be separated into two parts. Since we are dealing with a time-dependent 
system, we must specify the time for which a particular separation holds. 
We use the following definition: 

The stable part of phase space consists of the points {xi,X2), for which 
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the orbit 7(xi,X2,0) starting in (xi,X2) at t = remains bounded for t 
going to infinity. All other points define the unstable part of phase space. 

Clearly, a point {xi,X2) can only be contained in the stable region if it 
lies within an 0(e) neighbourhood of the area bounded by the homoclinic 
orbit of system (^). If this is not the case, 7(xi,X2,0) will reach the upper 
branch of the unstable manifold of the saddle point of system (P) in a finite 
time and clear off to infinity. We must not overlook the orbits starting close 
to the lower branch of the stable manifold of the saddle point of system (^), 
which can reach the just described 0(e) neighbourhood too. It will turn out 
that although this region may look small, it produces the major part of the 
stable region. 

These considerations help us locating the boundary of the stable region 
approximately. The location of the boundary of the stable region separates 
the part of phase space in which all orbits diverge to infinity ((xi,X2) — > 
(+00, +00)) from the part of phase space in which the orbits tend to circle 
around the origin for t going to infinity, so if we expect to see effects of 
the vanishing of the asymmetric potential somewhere, it is just within this 
boundary. 

The key step in analyzing system (^) is performing the transformation 



yi = a{et)xi (6) 

The idea behind this transformation is to try to fix the "slowly moving 
saddle point" of system @. Demanding that yi = y2, we arrive at the 
system 



= y2 

- yi + + 2e 2/2 + e ^ ^iet) ^^{^ j ^1 

where a'(et) stands for ^^^^ ^ ^ and similarly for a" {et). 

By transformation ^ the slow time-dependence has moved to 0(e) 
terms; still, system ^ looks more complicated than system However, 
we will be able to neglect the O(e^) term in most of our calculations. We 
should also note that system (0) is not Hamiltonian anymore, since we have 
applied a non-canonical transformation. Indeed the 0(e) term is a friction 
term, causing the origin (7/1,7/2) = (0)0) to become an attracting focus in- 
stead of a center. In the analysis of system (|^) we start with a special choice 
of a(et). 



2.1 The special case a{et) = e 



-tt 



We will first calculate the location of the boundary of the stable region for 
the special, but physically important case 
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a{et) = e"^* 



(8) 



We will show later that the general case does not differ much from this 
special case. With the choice (|8|) for a{et), system (0) reduces to 



y2 = -yi + yt- - e^yi 

It is remarkable that for this special yet interesting choice of a{et), our 
system becomes autonomous, which reduces the calculations because the 
dependence on the initial time has vanished into the transformation @. 
We also note that we have succeeded in fixing the saddle point: The saddle 
point of system (^) is located in (1 + e^,0). 

The saddle point not being located in (1,0) as we intended would in- 
troduce a lot of extra small terms in our calculations. To avoid these we 
map the saddle point onto (1, 0) by substituting — > (1 + e^)yj, i = 1, 2, to 
obtain 



2/1 = 2/2 

2/2 = -il + e'){yi-yi)-2ey2 



2u.. „,2^ o.„, (10) 



So we have reduced the calculation of the boundary of the stable region of 
system (^) to the calculation of the (time-independent) region of attraction 
of system (0). 

It is easily seen (figure ^) that the region of attraction of system ( [lO| ) is 
bounded by the stable manifold of the saddle point. 

It is well known that generally the stable manifold of a perturbed system 
(with parameter e) lies in an 0{e) neighbourhood of the stable manifold of 
the unperturbed system (with e = 0). The unperturbed system 



= , (11) 

= -2/1+2/1 ^ ' 

is simple and totally understood. It has a first integral E(e = 0) where 



E{e = 0) = ^yi + ^yj - ^y? (12) 

and the unstable manifold coincides with the homoclinic orbit E(e = 
_ 1 



0) 
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Figure 2: The stable and unstable manifold of system (|To| ) with e = 0.1 



Using E{€ = 0) in our calculations for the perturbed system would in- 
troduce some higher order terms. Instead, we extend the definition of E 
with suitable O(e^) terms which cancel these terms. Again, this is only for 
calculational convenience. 



E = lyi + l{l + e')yl-l{l + e')yl (13) 

It is instructive to combine figure |2| with the homoclinic orbit of the 
unperturbed system (11), which produces figure 

We will now approximate the location of the stable manifold of the sad- 
dle point of system (jl^) by calculating the variation of E along the stable 
manifold. Since this variation is an 0(e) effect, we may use the unperturbed 
stable manifold in this calculation, which involves elliptic functions. From 
this variation of E along the stable manifold, we can deduce the location of 
the perturbed stable manifold. 

If we follow the flow along the stable manifold from a point (yiO) 2/20) to 
a point (yii,y2i) we get: 



dE r /"S/ii 

— = -2eyl ^AE= -2eyldt = -2e / 7/2^2/1 (14) 

The integral appearing in this expression has to be calculated with 0(e) 
precision, which allows us to substitute the explicitly known orbits of the 
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Figure 3: The homoclinic orbit (represented by the thin hne) of system (|TT| ) 
added to figure ^ 

unperturbed system (p!l|). These orbits (1/2(2/1)) are readily obtained from 
the definition of the first integral ( |T2|) . 

To calculate the variation of E along the upper branch of the stable 
manifold, we take {yu, 2/21) = (1, 0) and get, after some analysis as indicated 
above: 

which is valid for — | < yi < 1 and y2 > 0. 

For yi < — ^ we have E = O(e^) and therefore we get: 

E{y^) = l + ^e + Oie') (16) 

which is valid for yi < — ^ 

Taking (yii,y2i) = (-^,^21) we get: 

6 \5 15 J ^^^^ 
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which is vahd for — ^ < yi < 1 and 1/2 < 0. The special form of the 
error term arises from the fact that the homoclinic orbit is only an 0(-y/e) 
approximation of the stable manifold for yi close to 1 and negative 2/2 (just 
under the saddle point). This follows from the analysis in Haberman and 
Ho i. 

Taking (yii,y2i) = (l,y2i) we get 



EM = . + 2. I . + — -^)P^. + 1>"^^' ) + 0(.^ log.) 




(18) 



which is valid for yi > 1 and y2 < 0. 

To calculate the variation of E along the lower branch of the stable 
manifold, we take (yii,y2i) = (1,0) and making use of the expressions for 
the explicitly known lower branch of the stable manifold of the unperturbed 



system (11) we find 
Bin) = 1(1 + a + 2. f ^ + v/3(..-2)(2,. + l)'»/n ^ _ 
which is valid for yi > 1 and y2 < 0. 

So, we have now calculated the variation of E all over the stable manifold 



of the saddle point of system (|lOD . What is left to do is to deduce the location 
of the stable manifold itself from this variation, which is not very hard. 
Given a value of yi, one first calculates the corresponding value of E 



using the appropriate formula given above. Using the definition of E (13), 
one calculates the corresponding value of y2. This amounts to solving a 
third order polynomial, which can even be done explicitly. 

In particular one can compute the intersection of the stable manifold 
with the yi-axis, which occurs (approximately) in (— ^ — f^'^)' 

2.2 The boundary of the stable region for arbitrary a(et) 

We now return to the discussion of the general system (0) . It turns out that 
the analysis is essentially the same as for the special case a{et) = e"*^* 

We claim that the behaviour of system is (with a certain error) 
described by the system 



yi = 2/2 

,2,0 (20) 

2/2 = -yi + yi + 2e— — -y2 
a(0) 

as far as the location of the boundary of the stable region is concerned. 
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The idea behind this statement is that if an orbit of system starts at 
a distance 0{5) inside the unstable region, it will reach the upper branch of 
the unstable manifold after an interval of time 0{log5), because it has to 
pass the saddle point at a distance 0{S) (sometimes twice). 
Using Gronwall's inequality, it is easy to show that the orbit of system (|20| ) 
starting at the same initial point, will diverge at most from the 

exact orbit. 

Since we know the boundary of the stable region with precision 0(e^ loge), 
we must take 6 to be larger than O(e^loge) for our calculations to make 
sense. 

Consequently, the orbit of system (pO|) will diverge at most o(l) from the 
exact orbit and will thus diverge to infinity too. 

Thus, a starting point (yi, 2/2) lying more than 0(e^ log e) inside the unstable 
region produces an orbit diverging to infinity both in system (^) and in 
system (|20|). 

We can apply exactly the same argument to the stable region, which 
proves that the boundary of the stable region of system (|^) coincides with 
the boundary of the stable region of system ( pO|) up to 0(e^ loge). 

So we have the important conclusion that, to calculate the boundary 
of the stable region of system (^), we can use the formulas derived for the 
special case a{et) = e"*^* with e replaced by — ^^e- 

3 Averaging inside the stable region for arbitrary 

a{€t) 

Knowing the location of the boundary of the stable region we proceed to 
study the stable region itself (the unstable region is clearly not very inter- 
esting). We have to do this study in two parts in which we consider the 
interesting dynamics which takes place close to (we will make this more pre- 
cise) the boundary of the stable region (i.e. in the boundary layer) and in 
the inner domain. At a safe distance from the boundary layer, system (0) 
will behave more and more like a harmonic oscillator. The natural way to 
approach such a problem is to apply the theory of averaging. 

3.1 Averaging in the inner domain 

Averaging in the vicinity of the origin is a simple exercise involving averaging 
over harmonic functions. This is not what we have in mind; we shall average 
over a part of the inner domain as large as possible. This involves averaging 
over elliptic functions. 
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3.1.1 Calculation of the averaged equation 

To perform averaging, we need one or more quantities with a small (0(e)) 
time derivative, i.e. which depend slowly on time. A natural candidate for 
this quantity is the exact integral (jl^) of the unperturbed system, for which 
we have 



dt a{et) \ \ a{et) a{et) J J 



(21) 



To be able to average this equation, we have to put restrictions on a{et): 



«(0 



is bounded for all positive ^ 

(22) 

is bounded for all positive ^ 



Most decaying functions of interest satisfy these restrictions. Functions 
decaying extremely rapidly, such as a(^) = exp(— exp(^)), do not satisfy 
these restrictions. But since a(^) decays very rapidly, we can safely put a(^) 
equal to zero for all ^ bigger than some for which a(^o) ^ 1) without 
affecting the dynamics. Other examples of functions which do not satisfy 
(I22I) are functions which vanish in a finite time like a(^) = 1 — ^. Again we 
can restrict the time span such that this poses no problem. 



To average equation (21), we consider r = et as an independent variable 
and add the equation 



f = e (23) 
Since we only have to average the 0(e) part of equation (pi]), we have 



to average 1/2 along a periodic orbit of the unperturbed system (|llD. This 
amounts to calculating the integral of 7/2 along the periodic orbit and involves 
the period of the periodic orbit. This is in the spirit of averaging as for 
instance presented in Sanders en Verhulst 

To calculate / y^'^t, we make use of yi = y2, which reduces the calcula- 
tion to the action / y2dyi. The functional dependence of y2 on yi for the 



unperturbed system can be retrieved from the exact integral (12) and is the 

square root of a third order polynomial. 

Using this, we find that we also need this standard integral 



rb I / 1 3 6 

j \l {x - a){h - x){c - x) = — \/67r(6 - a)Vc - 02-^1 (~2'2'^' c" 



24) 
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with 2-^1 the hypergeometric function, which holds when a <b < c. 
The a, b and c are the exacts roots of a third order polynomial and are thus 
awkward expressions even for our simple unperturbed problem. Surprisingly, 
the combinations b — a and c — a reduce to manageable expressions: 



/I TT 

b — a = VSsin ( - arcsin(12£' — 1) H 



c — a 



V^cos (- arcsin(12£' — 1) 



(25) 



Substituting all this we get 



/ 



y2dyi =2^\/67r3sin2 Q arcsin(12^ - 1) + ^ ) x 



VScos ^- arcsin(12i? — 1)^ x ^^g-j 

I 3 sin arcsin(12S - 1) + |^ 
X 2F1 I --, -,3, 



2' 2' ' cos arcsin(12S - 1)) 

The factor 2 arises because we have to integrate once from 6 to a and 
once from a to b. 

To calculate the period of the periodic orbit of the unperturbed system, 
we apply the standard technique of separation of variables to the exact 



integral (12). This leads us through a calculation similar to the one above. 



resulting in: 



period =2\/6 



\/3cos ^1 arcsin(12£' — 1 

xK I ^ 

cos ( g arcsin(12£' — 1) 



(27) 



where K is the complete elliptic integral of the first kind. 
We finally obtain the averaged equation by dividing equation (26) by 
equation (|^) and adding some extra factors from equation (^T|): 



a{€t) period 
a{et) ^ ' 



11 



It does not add much to the understanding of the problem to write down 
the averaged equation in it's full form. That is why we omit this. All that 
matters is that the right hand side is an explicitly known function A{E) of 
E, which we can approximate to arbitrary precision, and of time. 



3.1.2 Validity of the averaged equation 

Since the averaged equation is an approximation of the exact system (^), 
we have to address the question of the accuracy of this approximation, on 
which timescale it holds and where in the stable region. 

We expect that the closer we start to the homoclinic orbit of the unper- 
turbed system, the less accurate the averaged equation will become. The 
dynamics splits up in two qualitatively different time intervals, the first in 
which the orbit slowly separates from the homoclinic orbit and the second 
in which the orbit slowly spirals towards the attracting origin. 

We start with the first time-interval. As we will show in section 3.2 , 
apart from a sub-boundary layer of size 0(exp(— i)), this time-interval has 
a size of 0(|) independent of the initial distance from the homoclinic orbit. 
The total error introduced by the averaging process in the first time-interval 
is of O(eTo) (To is the period of the unperturbed orbit corresponding to 
£'(0), the initial value of £"). A short explanation of this estimate is given 
in section ^. 

For the second interval we can make use of the standard averaging the- 
orems, from which we get that the introduced error on the second interval 
is of 0(e) and that we are allowed to extend the second interval to infinity, 
because all orbits are attracted to the origin (see Sanders en Verhulst 0, 
chapter 4). 

This attracting property of the orbits also implies that the error introduced 
from the first does not blow up. Therefore, the total error introduced by 
the averaging process is of O(eTo) valid for all time. As we will also show 
in section 3^, for orbits starting close to the homoclinic orbit of the un- 
perturbed system £"(0) = we have that Tq is of order — log(| — E{0)), 
which implies that the averaged equation can be used to approximate the 
dynamics up to a distance of 0(exp(— i)) from the homoclinic orbit of the 
unperturbed system. More quantitative details on the boundary layer will 



be given in section 3.2. We will also see in section 3.2 that the averaged 



equation indeed breaks down when we approach the boundary layer. 



3.1.3 Analysis of the averaged equation 

We now turn to the analysis of the averaged equation (p8|). The first thing 
one should notice is that the effect of the decaying function a(et) can be 
removed from the equation by transforming to the new time r 
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I r=-llog(a(.*)) 

i a{€t) = e-'^ 

Note that this transformation reduces to the identity transformation in 
the special case a{et) = e~^*. 

It is remarkable that, given condition (p2[), it is not important at all how 
o(^) decays to zero, the dynamics of the system does not change, apart from 
a rescaling of the time axis. 

Applying this transformation produces the autonomous, 1-dimensional 
system 



dE 

— = -eA{E) (30) 

We can solve this system explicitly by separation of variables, but un- 
fortunately we do not have a primitive of in the form of an elementary 
function. 

But we can draw some important conclusions from this system, of which 
the most important one is the existence of an adiabatic invariant: As 
noted, it is always possible to solve system (|^), which gives the solution 
E = E{E{0), er) as a function of the initial condition and slow time. Again, 
in principle one can solve this equation for E{0) as a function of E and er. 
Inverting the time transformation (^), one finds E{0) as a function of E 
and et. Since £"(0) is obviously time- independent, we reach the conclusion 
that 

There exists a global adiabatic invariant inside the homoclinic orbit of 
the unperturbed system with the exclusion of an exponentially thin boundary 



layer, valid for all time, determined by equation (3t). 



For special cases we are able to produce these calculations explicitly, 
which we will show now. To understand these cases well, it is important to 
know how / y2dyi, the period and A depend on E. This is shown in figure 

I 

It is clear that / y2dyi depends almost linearly on E throughout the 
entire interval. This is understandable, since it is similar to the dependence 
of the area of a disk on its radius. What is not transparent is that the 
derivative of this function goes to infinity as E goes to | , but so slowly that 
its integral remains bounded. 

The period is close to 2tt for small E as it should be, because in this region 
the unperturbed system behaves nearly like a harmonic oscillator with uj = 
1. When E goes to g, the period goes to infinity, because the orbits are 
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Figure 4: The dependence of / y2dyi, the and A on E 

approaching the saddle point, in the neighbourhood of which they will stay 
a long time for each passage. 

The quotient of the two, A, shows the linear behaviour of / y2dyi for small 
E, because the period is almost constant. However, A has a maximum 
(0.248320 ...) at E = 0.152640 . . . , after which it rapidly drops to zero. We 
could have predicted that A is small for E close to ^, since all the time the 
orbits are close to the saddle point, the righthand side of equation (^) is 
small (2/2 ^ 1), resulting in a small average. 

3.1.4 The adiabatic invariant 

We now turn to the calculation of the adiabatic invariant for E{0) small (we 
will make this more precise later on). 

To approximate the adiabatic invariant, we perform a Taylor expansion of 
A{E) around 0. We note that the hyper geometric function forces us to use 
^/E as expansion variable instead of just E. However, it turns out that the 
coefficients in front of the non-integer powers of E are equal to zero, at least 
to fifth order. After a long calculation we arrive at the following expansion, 
valid for < ^ < i 



A{E) = 2E- ^E^ 



155 
'54 



E' 



61135 
5184 



E^ 



825409 
15552 



E^ + OiE^ 



(31) 



To approximate the adiabatic invariant, we truncate the series after the 
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second order terms, since we are interested in the first non-trivial deviation 
from a slowly attracting focus. Substituting this quadratic expression into 
the averaged equation (^) we get 



^ = -2eE + -eE^ (32) 
dr 6 



which is easy to solve giving 



Bi„) = , '-^ (33) 

(2 - |E(0)) e'" + |E(0) 

From this we readily obtain the adiabatic invariant 



2-1^(0) 2-p(er) 
We are now able to specify what we meant with £'(0) small. Since we 

— 2 

have neglected 0{{E{0)) ) terms, we have introduced a new error of order 

— 2 

(£^(0)) in the approximation of the solution. Since we do not want this error 
term to dominate the error introduced by the averaging process (0(e)), we 
take E{0) to be 0{e^/^). 

Expanding the adiabatic invariant around E = 0, we see that the first 
non-trivial correction to the slowly attracting focus (with adiabatic invariant 
£"(0) = Ee^'''^) is given by ^E'^e^'''^ resulting in a slightly slower collapse 
onto the origin (1/1,2/2) = (0,0). 

These arguments hold for the {yi , 2/2 ) phase space only. To extend them to 
the original {xi^X2) phase space, we have to invert the time-transformation 
(29) and the phase space transformation (^), after which we obtain the 



adiabatic invariant in the (xi,X2) phase space: 



2,a[et)xl - 2a{etfx\ + Qea'{et)xiX2 + 2>a{et)xl 
72a(et) - 15a(et)^xf + 10a(et)^xf - 30eo(et)^a'(et)xiX2 - lba{etf x^^^^^-^ 

We include this rather lengthy expression, because it reveals an impor- 
tant phenomenon: Due to the cross-terms X1X2, the level curves of the 
adiabatic invariant for a fixed time "resemble" ellipses, of which the long 
axis and the short axis differ by an 0(e) amount, and which are rotated 
around the origin, causing asymmetry with respect to the yi and y2 axis. 
We did expect this for finite time, but this behaviour persists when we let t 
tend to infinity. Put in other words, when t goes to infinity, our dynamical 
system (^ becomes symmetric (with respect to xi and 3:2), but the level 
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Figure 5: A few level curves of the adiabatic invariant for t fixed at infinity 

curves of the adiabatic invariant remain asymmetric. We have reached this 
important conclusion: 

The evolution of an ensemble of phase points towards a symmetric po- 
tential will show significant (i.e. 0(e)) traces of its asymmetrical past, for 
all time. 

So there is a sort of hysteresis effect present: although the system be- 
comes symmetric, it still "knows" that it was asymmetric in the past. 
We note that this phenomenon is not present in the (y 1,1/2) phase space, 
where the level curves of the adiabatic invariant are symmetric with respect 
to the yi-axis, but is introduced by the phase space transformation alone. 

To demonstrate this phenomenon visually, we have to take e not too 
small, so we took e = i. Fi gure ^ shows a few level curves of the adiabatic 
invariant for a(et) = e~^^ and t fixed at infinity. The asymmetric effect is 
clearly present. 

As explained before, we expected to see effects of the slowly decaying 
asymmetry in the neighbourhood of the boundary layer separating the stable 
and unstable region, but now it turns out that there are effects (0(e)) close 
to the origin too. 
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Figure 6: The structure of the boundary layer 
3.2 Approaching the boundary layer 

We study the approach to the boundary layer, which is an o(l) domain near 
the homoclinic orbit and the stable manifold. 

More precisely, the boundary layer can be divided into three regions (see 
figure ^). The first region consists of the phase points which are between 
0(exp(— -i-)) and o(l) inside the homoclinic orbit of the unperturbed system. 
It is in this region that the averaging technique slowly loses its validity, as 
explained in section ^. We will call this region the o(l) boundary layer. 

The second region consists of the phase points which are within an 
0(exp(— ■^)) neighbourhood of the boundary of the stable region. Orbits 
starting in these points will pass the saddle point (2/1,1/2) = (1; 0) on at least 
a J timescale (which requires special attention), after which they will enter 
the third region. We will call this region the 0(exp(— i)) boundary layer. 

The third region consists of the remaining phase points in the boundary 
layer, which is a strip with an 0(e) width. Orbits starting in this region 
will enter the first region on an 0(1) timescale, which allows us to use the 
unperturbed orbits inside this region. We will call this region the 0(e) 
boundary layer. 

The inner region, finally, consists of the phase points inside the stable 
region but outside the boundary layer. 

Using the same approach as in the previous subsection, we are able to 
study the adiabatic invariant everywhere in the inner region and the o(l) 
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boundary layer. The general idea is to expand the averaged equation around 
a certain value of E, in the neighbourhood of which we want the study the 
adiabatic invariant. This can be done to any desired precision. For low 
order expansions, it is possible to integrate the resulting equation explicitly. 
For high orders, one has to use numerical methods. 

Approaching the boundary layer, there are two more special values of 
E which we will study now, knowing the value of E corresponding to the 
maximum of figure (^) and the maximum value E = ^. 

The first special value can be computed numerically, giving Emax = 
0.1526396.... Expanding the averaged equation again to second order 
around Emax, we arrive at 

-cie + C2e{E - Emaxf (36) 
and C2 = 64.73966 . . . , which is easy to solve 



dE 
dr 

with ci = 0.2483204 . . 
giving 



E - Emax = y^tanh {y/c^{-eT + Ie^^J) (37) 

where /p is an integration constant determined by the initial con- 
dition. Solving lE^ax this equation, we arrive again at the adiabatic 

invariant: 



hma. = T^J^^^t^^^ (V - ^--)) + (38) 

Note that these equations hold only in an 0(e^/^) neighbourhood of 
Emax, and therefore on an e"-^/'^ timescale. For instance, it does not make 
sense to take the limit r oo. Although this limit does exist, its value is 
obviously wrong. Therefore the tanh should be regarded only as the first 
non-trivial correction to the linear time evolution of the adiabatic invariant 
around Emax- 

The second special value of E, g, is much more interesting and much 
more tricky, since it lies outside the domain of validity of the averaging 
process. We can however still expand the averaged equation around this 
value, because the part of the boundary layer inside the homoclinic orbit 
of the unperturbed system (0(exp(— ■^))) is much smaller than the domain 
of validity of the expansion 0(e^/^). Therefore, we are allowed to use the 
results of this expansion, but only outside the boundary layer. Indeed we 
will see that the results of this expansion inside the boundary layer are not 
correct. 
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Expanding the averaged equation around g is not simple, because the 
hypergeometric function has an infinite derivative at this point, and the 
elhptic integral (the period) is unbounded at this point. 

We break up this calculation by expanding equation (^6|) and equation 
(|27|) separately. After a straightforward calculation, we arrive at the follow- 
ing expansions: 



y2dyi =- - 72 I 1 - log 



E 



E 



+ O los 



72 
1 



72 



E\{--E) 



(39) 



period = — log 



72 



12 26 + 5 log 



72 



E 



72 



0\\og[\-E\ {\-e{ 



(40) 



Substituting these two expansions, we obtain for the averaged equation 
(with O (e(| — E)^^ terms neglected) 



d7 



12 



Slog 



72 



+ 144e 



V 



26 



log 



72 



5 log^ 



72 J 



72 



This equation is too complicated to be solved analytically. However, if 
we neglect the O ^e(g — E)^ term too, it is again possible to calculate the 
adiabatic invariant explicitly: 



h,„ = (i^) log (i^) - (i^) + (42) 

Note that this adiabatic invariant is only valid on an timescale, since 

\ — E will become 0{^/e) on this timescale, causing an extra error of 0(e). 
At this point we are able to make some important remarks: 

• Every orbit starting inside the o(l) boundary layer will collapse onto 
the attracting focus (yi,y2) = (0,0) on an ^ timescale, independent 
of the initial distance from the homoclinic orbit (collapsing onto the 
origin in the (yi, 1/2) plane is equivalent to circling around the origin in 
the (xi, X2) plane). This follows directly from the adiabatic invariant 
(|4^), which forces the orbits away from the boundary layer on an ^ 
timescale. 
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• The averaging process breaks down in the small strip between the o(l) 
boundary layer and the homoclinic orbit of the unperturbed system, 
like we expected it to. If the averaging process would be valid there 
too, every orbit starting there would collapse onto (j/i,?/2) = (0,0) on 
an i timescale, which would imply that all these orbits stay within 
a certain bounded neighbourhood of the origin in the {xi^X2) plane. 
This cannot be true of course, because an orbit starting very close to 
the saddle point {xi,X2) = (1,0) inside the homoclinic orbit, will end 
up arbitrary far away from the origin in the {xi,x-2) plane. 

• The leading order behaviour of the period near the homoclinic orbit 
is given by — log ^-^72"^ • This is the cause of the break-down of the 
averaging process, since averaging is only valid if the period is o(^). 



4 The boundary layer 

After the previous study of the major part of the stable region, we will turn 
our attention to the remaining part of the boundary layer. Since the o(l) 
boundary layer is covered by the previous section, we only have to study 
the 0(e) and 0(exp(— i)) boundary layers. As we explained in the previous 
section, we cannot use the theory of averaging for this study. 

We treat the 0(e) and 0(exp(— i)) boundary layers simultaneously. The 
only difference between them is that orbits starting inside the 0(exp(— i)) 
boundary layer will pass the the saddle point (2/1,^2) = (1)0) on ^it least a 
1- timescale, which results in an arbitrary large circle in the {xi,X2) plane 
as t tends to infinity. However this does not require a separate treatment. 

It is important to note that orbits starting inside the 0(e) boundary 
layer will remain within an 0(1) neighbourhood of the origin in the {xi,X2) 
plane. So, although the 0(e) boundary layer appears to be larger than the 
0(exp(— ■!)) boundary layer, it is in fact much smaller, because the latter 
has to fill up the rest of the (xi, X2) phase space. 

To study the boundary layer, we can use the same method we used to 
compute the position of the boundary of the stable region, because the orbits 
in the boundary layer are close to the homoclinic orbit of the unperturbed 
problem. So, to calculate the orbits to 0(e) precision, we are allowed to 
substitute expressions for the homoclinic orbit into the 0(e) contributions 
to the dynamics. 

This way we get again a two stage process. The first stage is governed by 
the homoclinic orbit of the unperturbed system. After the orbit has entered 
the domain of validity of the averaging process, the orbit collapses onto the 
origin on a ^ timescale. 

Since the existence of an adiabatic invariant was of great help in our 
study of the inner region, we prefer to extend that approach to the boundary 



20 



layer. We expect an adiabatic invariant to be present in the boundary 
layer too, because we are studying a Hamiltonian system which depends 
adiabatically on time. Finding this adiabatic invariant is generally very 
hard in regions where the unperturbed system has non-periodic solutions 
(in our case, outside the homoclinic orbit). 

The straightforward way to find the adiabatic invariant is to perturb 
the energy of the unperturbed system in such a way that its time-derivative 
becomes O(e^). So we are looking for an adiabatic invariant of the form 



Iu{yi,y2,et) = E{€ = Q)+eg{E{€ = Q),yi,et) (43) 
where E[e = 0) is given by ([l^). 

By demanding that the time-derivative of has a zero 0(e) contribu- 
tion, one normally arrives at a first order linear PDE for the function g. 
With a little bit of foresight, we choose the first argument to be orthogonal 
to the characteristic lines of the PDE, which is why we arrive at a first order 
linear ODE for the function g. 

By using Gradshteyn and Ryzhik |5| intensively we derived this explicit 
expression for the function g: 



/^/ N N 8 [2a'{et) 



{ (-!-{ + ^e) \/(«' - I'f + + (»' + f-c'jA - Ik] 



a = arccos 



Note that E(a, r) is the elliptic integral of the second kind and not the 
energy. The real numbers a, b and c are related to the roots of a third order 
polynomial in this way 



2E{e = 0) - yj + = |(2/i - a)((yi - a - bf + c^) Vyi 

(45) 

So a, b and c are functions of E{e = 0), with a < 0, 6 > and c > 0. We 
want to make four remarks with respect to this formula: 
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• We have derived the adiabatic invariant outside the homochnic orbit 
(but inside the stable region) of the unperturbed system. It is however 
not possible to calculate the adiabatic invariant in the inner region us- 
ing the same procedure, since the characteristic lines are closed curves 
in the inner region, which prohibits the PDE to have a solution. In- 
deed, the adiabatic invariant we found previously for the inner region 
has an 0(e) time-derivative. 

• determines the dynamics inside the boundary layer completely. 
This follows easily from d{Ibi) = 0. 

• Ibi is symmetric in y2. Transforming back to the {xi,X2) plane intro- 
duces again the cross-terms xiX2 in the adiabatic invariant which do 
not vanish for t going to infinity. 

• We now have an adiabatic invariant throughout the entire stable re- 
gion, with the exception of the very thin (0(exp(— i))) region between 
the homoclinic orbit of the unperturbed system and the o(l) boundary 
layer. This is not a problem, since we can approximate the dynam- 
ics inside this strip with transversal orbits which introduces only an 
0(^ exp(— i)) error. This trick solves the problem of matching the 
two adiabatic invariants at the same time. 

We would like to visualize the dynamics going on inside the boundary 
layer. Density functions are not very useful for this, since our system is 
Hamiltonian which implies area conservation. This is well known for au- 
tonomous and time-periodic Hamiltonian systems. To prove it for general 
time-dependent systems, one introduces a new independent variable equal 
to the time (making the system autonomous). Applying Liouville's theorem 
proves the desired result. 

Note that the conservation of area implies that the area of the "tongue" 
of the boundary layer is infinite, since it has to fill up the entire {xi,X2) 
phase space in the end. So our "thin" boundary layer is in fact the largest 
part of the stable region. 

To study the dynamics inside the boundary layer, we therefore choose 
to look at the evolution of the rectangular box around the homoclinic orbit 
of the unperturbed system, as depicted in figure |^. 

Since all orbits starting inside the box will remain inside the (evoluted) 
box, we only have to study the boundary of the box. Moreover, we only 
have to study those points of the boundary lying inside the stable region, 
since all other points clear off to infinity on an 0(1) timescale. Therefore 
we only have to study the bottom boundary (b) of the box. 

So by studying only a very limited set of phase space, we will gain infor- 
mation about all orbits starting inside the box, i.e. both inside the domain 
of validity of averaging and inside the boundary layer. 
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Figure 7: The rectangular box around the homochnic orbit of the unper- 
turbed system. 

For numerical reasons, we followed the evolution of the bottom bound- 
ary of a different (but similar) box in the {xi,X2) phase space, namely the 
straight line between (xi,X2) = (0,-2.5) and (xi,X2) = (5,-2.5). The nu- 
merical results are shown in figure |^. We took e = 0.1 which forced us to 
take steps along the boundary as small as 10~^^ to generate the last sub- 
figure. This is due to the fact that the most interesting dynamics takes place 
in an 0(exp(— i)) neighbourhood of the boundary of the stable region. 

The open area around the origin is the domain of validity of averaging. 
This is the part of phase space where the level curves of the adiabatic in- 
variant (figure ^ live. It is also clear to see the instantaneous saddle point 
moving from (1,0) to infinity. 

The points connecting the instantaneous saddle point with the domain of 
validity of averaging have started very close (0(exp(— i))) to the boundary 
of the stable region, passed the saddle point during a time-interval of O(^), 
after which they entered the domain of validity of averaging (in the (2/1,2/2) 
phase space). 

The effect in the (xi,X2) phase space is that the orbits end up circling 
around the origin outside the part where the level curves of the adiabatic 
invariant (figure H) live. The closer an orbit starts to the boundary of the 
stable region, the larger the radius of the circle it describes in the end. 

The effect of the area conservation is also nicely visible. Since the starting 
box has a finite area, the area inside the spiral must be finite too, which 
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Figure 8: The evolution of a straight line (part of the bottom boundary b) 
crossing the boundary of the stable region. 
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makes the spiral very thin. Note that from f = 4.7 on the curve going to 
infinity actuahy consists of two very close curves. 

Note also that the remaining (major) part of phase space has to be filled by 
the tail of the "small" tongue of the boundary layer which lies outside the 
box. 

5 Concluding remarks 

It is surprising that it is possible to give a fairly complete treatment of system 
which describes the evolution of a simple system with an asymmetric 
potential to a symmetric potential. The most remarkable result is that 
in the evolution towards symmetry as time tends to infinity, traces of the 
asymmetric past can be recognized in the solutions. 

System (^) is just a metaphor for simplified models with two degrees of 
freedom which exhibit evolution from asymmetry towards symmetry. In a 
forthcoming paper we shall discuss such higher dimensional problems using 
basically the same methods. 

In the discussion of the validity of the averaged equation, we have as- 
sumed that the reader is familiar with the proof of the standard averaging 
theorems (see for instance Sanders and Verhulst ||^]). In particular, we have 
used the straightforward extension of those theorems to periods depending 
on e: by rescaling the time- variable it is easily shown that averaging pro- 
duces 0(eT(e)) approximations, valid on a - timescale, as long as eT(e) is 
o(l). 
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